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The ground state phase of spin-1/2 J1-J2 antiferromagnetic Heisenberg model on square lattice 
in the maximally frustrated regime (J2 ~ 0.5Ji) has been debated for decades. Here we study this 
model by using a recently proposed novel numerical method - the cluster update algorithm for tensor 
product states (TPSs). The ground state energies at finite sizes and in the thermodynamic limit 
(with finite size scaling) are in good agreement with the state of art exact diagonalization study, and 
the energy differences between these two studies are of the order of 10~^ Ji per site. At the largest 
bond dimension available Z) (D = 9), we find a paramagnetic ground state without any valence 
bond solid order in the thermodynamic limit in the range of 0.5 5S J2/J1 ^ 0.6, which implies the 
emergence of a spin-liquid phase. Furthermore, we investigate the topologically ordered nature of 
such a spin-liquid phase by measuring a nonzero topological entanglement entropy. 



PACS numbers: 

Introduction 

The spin 1/2 J1-J2 antiferromagnetic Heisenberg 
model on a square lattice has drawn great attention 
for the last two decades owing to its close relation to 
the disappearance of the antiferromagnetic (AF) long 
range order (LRO) in the high- Tc superconducting ma- 
terials [U [2] , and has been proposed as a possible simple 
model to realize topologically ordered chiral spin-liquid 
state ini H] or spin- liquid state [SHE] . The Hamiltonian 
of this model is given by: 

i7 = Ji^Si-Sj + JaJ^Si-Sj, (Ji,J2>0), (1) 

where represents the nearest neighbor (NN) pair 

and represents the next nearest neighbor (NNN) 

pair. For convenience, we set Ji = 1 throughout the 
paper. It has long been believed that the frustration 
from NNN interaction competes with the NN one and 
drives the system through a quantum phase transition 
from an AF LRO phase to a magnetically disordered 
phase. In two extreme cases, the ground state phases 
of the model are well established: at very small J2, the 
ground state has AF LRO; and at very large J2, the 
system falls into two weakly coupled sets with magnetic 
susceptibility peaks at momentum (tt, 0) or (0,7r). In 
the intermediate coupling regime, quantum fluctuation 
is meant to destroy the AF LRO before the maximally 
frustrated point J2 = 0.5 of the classical model and es- 
tablishes a new paramagnetic phase. The nature of such 
a quantum phase is of great interest. 

Numerous efforts have been made using many dif- 
ferent approaches, such as the exact diagonalization 
(ED) ISHHli spin- wave theory [TSl [11], series expan- 



sion [TTI HH] , large-N expansion [19] , the coupled cluster 
method (COM) 120], variational methods (including short 
range resonating valence bond (SRVB) method) [2m23j . 
and the fixed- node quantum monte carlo (QMC) |24j . 
Especially, a series expansion calculation of a general 
magnetic susceptibility over different perturbation fields 
suggests that within the Ginzburg-Landau paradigm the 
type of phase transition from the Neel to paramagnetic 
phase is of first order [18j. However, the same general 
magnetic susceptibility calculated with a coupled clus- 
ter method suggests a second order phase transition '5D] . 
The debates about the phase near J2 — 0.5 are more 
or less the same. A fixed-node QMC study indicates a 
plaquette valence bond solid (PVBS) state whereas 
the series expansion argues for a columnar valence bond 
solid (CVBS) state |TB|. A relatively direct investiga- 
tion of the nature of the ground state order is using the 
SRVB approximation [22', where with another term J3 
included in the Hamiltonian, a PVBS state along the line 
of J2 + J3 — 0.5 is found. It is worth to mention that the 
recent TPS studies have resolved a CVBS state [IS], but 
with rather small bond dimension. 

In this paper, we revisit this problem with a TPS |26| 
ansatz for the ground state wave function, accessed by 
the recently proposed cluster update algorithm |27| , and 
reveal the answer to both questions. Up to D = 9, we 
observe a continuous phase transition from the Neel to 
paramagnetic phase at J2 = 0.5. We further investigate 
the nature of the paramagnetic phase by measuring var- 
ious VBS order parameters. Through finite size scaling 
(ESS), all considered VBS orders (including CVBS, stag- 
gered VBS (SVBS) and PVBS) scale to zero, implying 
the possible emergence of spin-liquid state. The expo- 
nentially decaying spin-spin and dimer-dimer correlation 
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FIG. 1: (a) The finite size {L = 16) ground state ener- 
gies at different bond dimensions D — 3,4, 5, 6, 7, 8, 9, and 
their extrapolation (from ground energies at bond dimensions 
D = 6, 7, 8, 9) to D — ^ oo Umit. (b) The finite size scaling of 
E{J2,L) obtained from (a) to the thermodynamic limit. 



functions found in the paramagnetic phase indicate the 
non-zero spin and magnon excitation gaps. Moreover, we 
further verify the topologically ordered nature of such a 
spin-hquid phase by computing its constant part of Renyi 
entanglement entropy. Comparing to the previous stud- 
ies in other systems [5S1 US] , our results provide stronger 
evidence for the existence of spin-liquid states in the spin 
1/2 J1-J2 antiferromagnetic Heisenberg model on square 
lattice. 



Results 

We divide the square lattice into four sublattices 
A,B,C,D and associate each sublattice with a different 
tensor. Such a choice of tensor product state ansatz aims 
at describing all possible VBS orders and studying their 
competing effects. We use the cluster update imaginary 
time evolution method to evolve from a random TPS to 
the ground state of the above Hamiltonian. Upon ob- 
taining the ground state TPS with bond dimension D 
up to 9, we evaluate the ground state energy, magneti- 
zation, dimer order and plaquette order parameters of a 
finite lattice with periodic boundary condition (PBC) of 
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FIG. 2: (a) Finite size energies E{J2,L) obtained as in 
Fig.[T]^a) and their thermodynamic limits E{J2) from the FSS 
in Fig. [ijb). The purple and blue dashed lines are the ED 
energies for lattice sizes L = 4, 6 respectively, (b) Energy ab- 
solute error per site for system size L = 6 as a function of 
J2/ J\ compared between the TPSs study at bond dimension 
D = 9 (red squares) and a VMC study (green dots) with a 
RVB wavefunction ansatz [23] • The blue triangles are energy 
differences between E{J2) from the TPSs study and those 
extrapolated from an ED study up to 40 spins |14| . 



size i X L for L = 4, 6, 8, 12, 16, 32 based on the renor- 
malization concept |30H32j . using the monte carlo sam- 
pling technique [32]. We note that the TPS obtained 
from cluster imaginary time evolution method is indeed 
a variational state of an infinite system, however, it is 
very difficult [33] to evaluate the physical quantities for 
large systems with PBC. Here we evaluate them on rel- 
atively small size systems with high precision and then 
extrapolate to the thermodynamic limit. 

We present the ground state energy for a finite lattice 
obtained from the TPSs of various bond dimensions D in 
Fig.jlja). A clear decreasing of the ground state energy is 
observed when increasing bond dimension D, especially 
at J2 G (0.45,0.6). We thus extrapolate the finite size 
energy to the _D — )• 00 limit with a formula 



E(J2,L,D)^E{J2,L) + c.jjD\ 



(2) 



here E{J2, L) and cj^ are the fitting parameters. A finite 
size scaling formula |14[ 134) is then employed to extrapo- 
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FIG. 3: (a) The FSS of the magnetization square to their 
thermodynamic hmits. The maximal system size is 32 x 32 
in this fitting, (b) The finite size magnetization M(J2,L) 
and their thermodynamic hmits M(J2) obtained from (a) as 
a function of J2 / Ji . 



late the ground state energy in the thermodynamic limit, 

E{J2,L)^E{J2)+djjL\ (3) 

here E{J2) and dj^ are the fitting parameters. The FSS 
results arc presented in Fig. [IJb). The extrapolated fi- 
nite size ground state energies together with their ther- 
modynamic limits are plotted in Fig. [2ja). The differ- 
ence between E{ J2) and the thermodynamics limits ob- 
tained from a recent ED study [14 is plotted in Fig.[2][b). 
A good agreement (within the order of 5 x per 
site) is reached between these two methods, whereas the 
rather large differences at J2 = 0.5, 0.55 are caused by 
the unreliability of the ED extrapolation at these cou- 
pling strengths [13]. Fig. [2jb) also compares the abso- 
lute energy error per site for system size 6x6 obtained 
from the TPSs study at bond dimension D ~9 with the 
RVB variational MC results f53| , a better achievement is 
obtained through the TPSs approach. 

In the MC sampling scheme, the staggered magnetiza- 
tion square is evaluated as 



1 ^ 



7V2 



(4) 



where 



Uj and {xi,yi) are the lat- 
tice coordinates of site i. We present the magnetization 



results for finite size lattices in Fig. |3] Here the data pre- 
sented are the results with the largest bond dimension 
available {D = 9) and without extrapolation over D, be- 
cause the magnetization is not a monotonic function of 
D. A FSS formula [T3J[33] is again employed to extrap- 
olate the magnetization square in the thermodynamics 
limit, 



M^J2,L) = M^{J2) 



(5) 



here J2), Cj^ and Cj^ are the fitting parameters. The 
FSS of the magnetization square is plotted in Fig. [sf^a) 
and the magnetization at finite sizes and their thermo- 
dynamic limits are presented in Fig. [3](b). We observe 
that the magnetization decreases to zero at « 0.5. 
This is clearly a sign of a continuous phase transition 
between the AF phase and the paramagnetic phase. To 
encounter a different asymptotic behavior in the param- 
agnetic phase, we employ a slightly different fitting for- 
mula for the finite size magnetization square J2, L) = 
]VP{Ji)+e^j^/ L^^^ . The power law (exponentially) decay- 
ing spin-spin correlation function at J2 = 0.5 ( J2 = 0.55) 
is another proof for the disappearing of the staggered 
magnetization on and after the critical point « 0.5 
(see supplemental material Fig.js]). 

To further determine the phase at the paramagnetic 
region J2 € [0.5,0.6], several order parameters are inves- 
tigated. We define the CVBS orders as 



ML 



Mi 



iV2 



N 



1 



iV2 



and the SVBS orders as 

1 ^ 



Mi 



N 



M^ = 



(6) 
(7) 

(8) 
(9) 



here ZJf = ^{x^.vi) ' S(2,,,y,)+s, (" = ^,2/) and 

Xi — Xj +yi~ Uj . We also define the P VBS orders as [! 



M2 = 



1 ^ 

^(_1)0..((Q^Q^.)_(Q^)(Q^.)), (10) 



Q JS12 

here Qi = ^{POi is the permutation operator that 

permutes the 4 spins on a plaquette. 

Fig. |4] and Fig. [5] demonstrate the CVBS, SVBS and 
PVBS orders respectively as a function of inverse lattice 
size L evaluated using TPSs at bond dimension D = 9. 
Through FSS, we find all of them scale to zero in the 
thermodynamic limit. This result rules out the VBS as a 
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FIG. 4: The CVBS order parameter (a) and (b) M^y as 
a function of inverse system length 1 /L fitted using a second- 
order polynomials, evaluated from TPSs of bond dimension 
_D = 9 at L = 4,6,8,16,32. 

candidate for the ground state of the paramagnetic phase, 
and the most possible state of the paramagnetic phase is 
a spin-liquid state. We further compute the real space 
spin-spin and dimer-dimer correlation functions (see sup- 
plementary material) and find all of them decay exponen- 
tially with distance in the paramagnetic phase, implying 
the gaped spin-liquid nature of the observed paramag- 
netic phase. Moreover, a metastable state (see Fig. [T] 
which has slightly higher variational energy) with coex- 
isting CVBS and SVBS orders is observed at D < A (see 
supplementary material), which is consistent with pre- 
vious study pS] and implies the emergence of the spin- 
liquid phase as a consequence of melting various VBS 
orders. 

To determine whether the observed paramagnetic 
phase is topologically ordered or not and what kind of 
topological order it is, we further compute the Renyi en- 
tanglement entropy S2{L) between two L x L subsys- 
tems by cutting a 2L x L torus into equal size halves. 
The entanglement entropy for a cylindrical geometry is 
presented in Fig. |6][a) . Through extrapolation we could 
isolate the the constant term of the Renyi entanglement 
entropy S2 of the paramagnetic phase. It has been shown 
that the translational invariant TPS (with a 2 by 2 unit 
cell) ansatz belongs to a special topological sector of the 
degenerate ground states and the Renyi entanglement en- 
tropy S2 will scale as S2 = aL — 7 where 7 is 
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FIG. 5: The SVBS order parameter (a) M|^, (b) M% and (c) 
the PVBS order parameter Mq as a function of inverse system 
length 1/L fitted using a third-order polynomials, evaluated 
from TPSs of bond dimension D = 9 at L = 4, 6, 8, 16, 32. 



a universal constant - the topological entanglement en- 
tropy. We find that in the paramagnetic phase 7 reaches 
a conspicuously large value between 0.8 and 1.3. One 
possible candidate for the paramagnetic phase, the gaped 
Z2 spin-liquid jH |40] , has a topological entanglement en- 
tropy equal to In 2 ~ 0.69 [35l |36], which happens to 
be close to the calculated intersection. As an additional 
check, we also calculate the Renyi entanglement entropy 
between a square of size Lx L and the rest of the system 
on a 16 X 16 lattice for L S [3:8]. The Renyi entan- 
glement entropy for such a disk geometry is presented in 
Fig. [6]Jb). We also find an nonzero intersection in this 
case, however, due to the non-universal corner entropy 
contributions, the actual value has a small difference. 
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FIG. 6: (a) The Renyi entropy S2(L) between two L x L 
subsystems by cutting a 2L x L torus into halves, evaluated 
from TPSs of bond dimension D = 9 at L = 4, 8, 16. (b) 
The Renyi entropy S2{L) between a L x L disk subsystem 
and the rest of a 16 x 16 torus, evaluated from TPSs of bond 
dimension D = 9 for L € [3:8]. 



the paramagnetic phase, and found strong evidence to 
support the possibility of being a spin-liquid. In princi- 
ple, the complete topological information characterized 
by the fractional statistics of quasiparticles can be read 
out through a renormalization scheme [2] for the TPS 
ansatz. A detailed study along this direction will be pre- 
sented in our future work. 



Method 

The cluster update imaginary time evolution 
algorithm 

The following is an illustration of how to construct the 
evolution operators for this Hamiltonian. We expand the 
evolution operator O ^ exp{— eJi(Si • S2 -|- S2 • S3) — 
2eJ2Si • S3} on 3 sites from the Trotter decomposition of 
the partition function. By writing exp(— eJS; • Sj) as 

Yl [cosh (eJ/4) 1, (g) Ij ~ sinh (e J/4) < (g) af] , (II) 

a 

where a = x,y,z, a" are Pauli matrices, and omitting 
higher orders of 0{e), one obtains 

6 = li I2 «> I3 - X! tanh(eJi/4)cr" ctj ® ^3 

a 

-^tanh(eJi/4)li (g) cr^ (g) 

a 

-^tanh(eJ2/2X(gl2«)CT^ (12) 



Conclusions and discussions 

In conclusion, we applied the cluster update algorithm 
for TPSs to study the frustrated spin 1/2 J1-J2 antiferro- 
magnetic Heisenberg model on square lattice. Limited to 
a cluster size 2 x 2, a rather large bond dimension D = 9 
is feasible. Through a finite D scaling and finite size scal- 
ing (FSS), our ground state energies at finite sizes and 
in the thermodynamic limit are in good agreement with 
the results from a state of art ED study J4j, with a dif- 
ference at an order of 10~'^ Ji per site. Through FSS, the 
staggered magnetization diminishes to zero at « 0.5, 
suggesting a continuous quantum phase transition at a 
critical point w 0.5. We further determined the na- 
ture of the paramagnetic phase using 3 sets of order pa- 
rameters: the columnar dimer order AfJ^ and M^y, the 
staggered dimer order M^^ and M^y, and the plaque- 
tte order Mq. At a large bond dimension D {D — 9), 
we found that all these order parameters scale to zero 
in the thermodynamic limit, which ruled out the colum- 
nar, staggered and plaquette dimer orders. We computed 
the constant part of the Renyi entanglement entropy for 



The above terms can be expressed as a matrix product 
operator (MPO) [12], 

3 

6= KB,,v,3)X,,0X,,®X,3 

Xo = i, Xl = a^ X2 = ^T^ x3 = fT^ 

vo = |0), 

Vj = a|i), (i== 1,2,3), 

Bo = |0)(0|+6|l)(l|+5|2)(2|+6|3)(3|, 

B, = c|0)(z|+c|z)(0|, (*- 1,2,3), (13) 

where are the vectors of length 4, are 4x4 matri- 
ces, Xi are operators acting on the physical index, and 
a, b, c are scalar variables. In order to correctly match 



the coefficients in front of each term in Eq. (12), a,b,c 
have to be chosen to satisfy ac = — tanh(e Ji/4), a^b — 
— tanh(eJ2/2), and |a|,|5|,|c| ^ 1. Thus the evolution 
operators on these sites are written as Oi = J2i vf (gX^, 
O2 = I]i g) Xi and 63 = ^ . Vi g) X^ respectively. 

We present the diagrammatic representation of the 
evolution operators Oi, O2 and O3 acting on sites A, 
B and C in a 2 x 2 cluster in Fig. [Tj^b) . The correspond- 
ing simple update scheme is sketched in Fig. [7][a). In 
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for passing their variational monte carlo data for com- 
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ERC grant QUERG, and the FWF SFB grants FoQuS 
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Note added. — After the completion of our work, we 
learned that Hong-Chen Jiang, Hong Yao and Leon Ba- 
lents obtained similar results [44j from high precision 
DMRG calculation. 



FIG. 7: (a) The simple update scheme, (b) The cluster up- 
date scheme with a cluster size 2x2. (c) A illustration of the 
MC sampling of the swap operator to measure the Renyi en- 
tanglement entropy S2 of a subset /3 for a 1 dimension system 
using 2 copies (represented by red and blue dots respectively) 
of the system. Dashed lines indicate summing over the spin 
degrees of freedom. W{ai,/3i) is the weight of the spin con- 
figuration {at, Pi) calculated using a tensor RG contraction 
scheme in 2 dimension. 



both cases, the complexity scales as D^, and there is no 
cumulative error. 



Monte Carlo sampling of Renyi entanglement 
entropy 

The Renyi entanglement entropy S2{L) can be mea- 
sured via Monte Carlo sampling of the swap operator 
using two copies of the system [53], as demonstrated in 
Fig. He), 



{Swap) 



W\a^)W\o2)Swap{a^,a2) 



Swap{ai,a2) 



W{ai,(32)W{a2,f3i) 
iy(ai,/3i)VF(a2,/32)' 



(14) 



where ai = (a.i,j3i), ai,/3i are the spin configura- 
tions of the left and right halves of the copy i, and 
W{ai,l3j) is the coefficient of the spin configuration 
{ai,/3j), which can be evaluated using a tensor RG con- 
traction scheme [3TJ [32] ■ 
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SupplementEtry material 



The spin-spin and dimer-dimer correlation functions 



We further calculate the spin-spin and dimer-dimer correlation functions at both critical point and deep in the 
paramagnetic region. Define 

C(r,0) = (S(o,o)-S(,,o)), (15) 

Cd,(r,0) = (£>^D^), (16) 

we found the spin-spin correlation (— l)'"C(r, 0) decays as a power law r~"= at the critical point with an exponent 
as = 1.31(3), however decays exponentially e"''/'^^ in the paramagnetic region with a spin correlation length = 
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2.6(1); and the dimer-dimer correlation {—lyCdxir^O) decays as a power law r~"'^ at the critical point with an 
exponent ad = 3.0(3), while decays exponentially e~'^l^'^ in the paramagnetic region with a dimer correlation length 
— 1.49(1), as demonstrated in Fig. [s] and Fig. [oj This indicates that the paramagnetic phase has both spin S — 
gap and spin 5 = 1 gap. 
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FIG. 8: The spin-spin correlation as a function of separation r at (a) J2 = 0.5 and (b) J2 = 0.55. At the critical point J2 — 0.5, 
(— l)'^C(r, 0) decays as power law r~°= with ~ 1.31(3), while at J2 = 0.55, (— l)' C(r, 0) decays exponentially as e^^^^" with 
a spin correlation length — 2.6(1). 



Observation of VBS order at small bond dimension D < 5 

In addition to the squared VBS orders, we also define the average of the dimer orders to detect if the symmetry is 
broken or not; 

Dd-x = Y,i~ir{Dt), Ddy = Y.{-ir{D\)- 

i i 



Fig. 10 shows the absolute value of the average order parameters as a function of J2/J1 at bond dimension D = A. 
We find that the paramagnetic ground state has nonzero Ddx,Ddy and Dgy of a magnitude 10~^, which indicates a 
mixed columnar and staggered VBS orders. However, after D > 5, all of them disappear when ground state energy 
further decreases (see in Fig. [l] for a comparison of energies at different bond dimensions D), implying the local 
minimal effects of these ordered states. The observed VBS orders for small D are a consequence of spontaneous 
symmetry breaking. However if increasing the bond dimension D, all average order parameters become zero, which 
indicates a restoration of certain broken symmetries caused by the restriction of bond dimension D. 
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FIG. 9: The dimer-dimer correlation as a function of separation r at (a) J2 = 0.5 and (b) J2 = 0.55. At the critical point 
J2 = 0.5, {—ly Cdx(r,Qt) decays as power law r"""* with ad = 3.0(3), while at J2 = 0.55, (— l)'^Cd2,(r, 0) decays exponentially 
as e"'"''*'' with a dimer correlation length = 1.49(1). 
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FIG. 10: The average columnar dimer order Ddx, Ddy and staggered dimer order Dsx, Dsy (defined in Eq. 17 1 as a function of 
Jij J\ at bond dimension D = 4 for system size L = 16. 



